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ABSTRACT 


The  theory  of  finite  Fourier  transforms  is  developed 
from  the  definitions  of  infinite  transforms  and  applied  to 
the  computation  of  convolutions,  correlations,  and  power 
spectra.  Detailed  procedures  for  these  computations  are 
given ,  including  listings  and  writeups  of  FORTRAN  subroutines 


introduction 


for  the  past  several  months,  E.  A  Flinn  t 
-a  1  have  been  examining  some  practical'and  '  "  ClMrt°Ut 

aspects  of  the  theory  of  Fourier  transf 

resulted  In  a  set  of  n  tra"^°™s.  These  efforts  have 

a  set  of  programs  for  performing 
series  based  on  the  Cool  „  Perf°rming  operations  on  time 

on  the  Cooley-Tukey  (References  1  51  v 
Fourier  franc-p  1/2;  hyper-rapid 

ier  transform  method.  usi  nr,  P 

seismic  array  data  such  as  the  LJL™'  °" 

correlations,  spectra,  and  digital  filters  h  C°nTOlUtl0nS ' 

by  fa-tors  of  three  or  four  and  s  •  S  ^  b0e"  Sf>eeded  UP 

of  this  report  is  to  c  eVen  ten'  Ths  Purpose 

eport  is  to  communicate  these  result.  • 

forward  manner  and  to  offer  some 

as  well  as  for  future  efforts  in  this  area  «•  ^  derivati°n 

listinqs  of  the  ”  a*  Writeups  and 

of  the  programs  discussed  here  are  included 
ces  to  this  report.  included  as  appendi- 

In  the  case  of  continuous  data  of  inflmt  , 

Fourier  transform  r,  -  lte  len?th/  the 

transform  pair  is.  usually  written  as* 


A  (uu)  =  — — 


2tt 


f  f(t 


-iout 


dt 


f(t) 


(1) 
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The  first  of  these,  going  from  time  to  frequency,  is  referred 
to  as  the  direct  transform  and  the  other  as  the  inverse  trans¬ 
form-  Sometimes  the  direct  transform  is  written  with  a  factor 
of  1  in  front  of  the  integral  and  the  inverse  with  a  factor  of 
1/2tt  .  These  are  of  course  equivalent  to  the  above  definition. 
Usually  the  quantities  of  interest,  such  as  spectra,  etc., 
involve  magnitudes  or  squares  of  one  transform  and  the  factor 
must  be  inserted  or  taken  out,  depending  on  which  definition 
is  used,  to  preserve  true  ground  motion. 

Two  drawbacks  of  these  definitions  for  digital  compu¬ 
tations  are  apparent:  First,  the  integrals  must  be  approximated 
by  sums  in  the  digital  computer,  which  implies  that  both  trans¬ 
forms  involve  sampled  variables.  Second,  the  infinite  limits 
on  the  sums  are  impossible.  Clearly  these  sums  must  truncated, 
as  they  do  not  in  general  converge  over  a  finite  interval. 

As  a  result  Fourier  transforms  as  such  are  never  really  com¬ 
puted  by  a  digital  computer.  Instead,  the  complex  samples  of 
a  direct  transform  »re  approximated  by  the  cosine  and  sine 
coefficients  of  Fourier  series  representation  of  the  input  data. 
The  definitions  for  these  are: 

00 

if  x(t)  =  ^  j^an  cos  (rrnt/T)  +  b^  sin  (rrnt/T)"|  ,  (2) 

n=0 

T 

then  a  =  ~  [  x(t)  dt  b  =  0  (3) 

o  T  J  o 

o 

T 

a  ~  ~  \  x(t)  cos  (nnt/T)  dt 
n  T  J 


2 


T 

2  r 

b  =  —  x  ( t )  sin  (rrnt/T)  dt 

n  T  J 

o 

If  N  samples  of  the  data  are  taken  at  equally  spaced 
intervals  At  =  T/N,  the  integrals  (3)  becomes  sums,  and  the 

f 

frequency  sum  in  (2)  goes  from  DC  to  the  folding  frequency, 
i.e.,  0  to  N/2T  .  The  equation^  are  then  written  ass 


N/2 

cos  (2it jk/N)  +  b^  sin  (2tt jk/N)  J 

k=0 

(4) 

N-l 

ao  =  N  I  x(j)  bo  =  ° 

j-0 


N-l  N-l 

ak  =  N  Z  cos  (2rr  jk/N)  bk  =  §  Z  sin(2njk/N), 

j-0  j=0 

where  t  has  been  replaced  by  jAt  »  By  now  defining: 

A(k)  ■  (a  -  i  b  )  ,  A (o )  =  a  (6) 

Z  K  K  O 

and  realizing  that  a  real  time  series  contains  only  real 
points,  (4)  can  be  written  as: 

N/2 

x(j)  =  l  A 00  exp  (2rri  jk/N)  .  (7) 

k=0 


A  great  deal  of  symmetry  between  the  two  transforms  can  be 


preserved  if  the  sum  in  O’  is  summed  up  to  N-l  .  Redundant 
points  in  the  spectrum  are  included  (since  the  transforms  are 
periodic)  but  the  computational  procedures  are  simplified 
It  is  also  convenient  to  split  the-  factor  of  1/N  appearing  in 
(5)  into  two  factors  of  1//n"  ,  one  in  front  of  each  transform 
By  defining  a  complex  number. 

w  *  exp(2ni/N)  ,  (8' 

the  two  transforms  can  now  be  written  as: 

f(j)  «-jk  O) 

a  no  wjk  .  do) 

It  can  be  shown  that  the  set  of  direct  Fourier  transform 
points,  between  DC  and  the  folding  frequency,  contains  the  same 
amount  of  information  as  the  real  data  series 5  The  transform 
includes  N/2  distinct  points,  which  with  the  DC  term  makes  a 
total  of  N/2  +  1  complex  points  Equation  (9)  shows  that  both 
the  DC  and  the  folding  frequency  point  are  purely  real?  thus, 
the  Fourier  transform  contains  (N/2-1 ) *2+2* 1  numbers.  This  is 
exactly  the  same  amount  of  information  contained  in  the  real 
time  series,.  It  also  suggests  that  the  existence  of  one  trans¬ 
form  should  imply  the  existence  of  the  other. 

If  there  are  N/2+1  non-redundant  points  in  the  direct 
transform,  then  the  sampling  interval  in  frequency  must  be 


N-l 

a  no  =  —  y 

j-o 

N-l 

I 


v'N 


k=0 
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(N/2T)/(N/2)  -  1/t  .  Thus,  the  product  of  the  time  and  frequency 
variables  is: 


•  T 

luut  =  i  2rrj  — 
J  N 


2rri 

N 


jk 


(11) 


This  equation  relates  the  arguments  in  the  two  exponentials, 
one  in  the  continuous  transform  and  the  other  in  the  finite 
transform  (Equations  1,  9,  and  10). 


3 •  TWO-AND  THREE-DIMENSIONAL  FOURIER  TRANSFORMS 

Two-and  three-dimensional  direct  Fourier  transforms  are 
seen  to  be 


N  -1  N  -1 


A(krk2)  = 


i_  *  f  ,  -j  k  - 

TT  l  l  X(  W  “j  1  1  “2 


-s/fo.N,,  ,  . 

12  J  =0  j  =0 
1  2 


^2k2  (12) 


and 


A  (kl  ,k2  ^3)  = 


V1  n2-i  n3-i 


- _  I  7  Y  x(j  ,j  j  )„ 

mm4-'  “-1  L‘  12  3  1 


123  31=°  32=°  j3=° 


W2  W3_ 3  3^3  (13) 


We  can  break  up  Equation  (12)  as  follows: 


n2-i 


a<VV 


=  l  B<kl'j2> 


(14) 
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This  calculation  requires  one-dimensional  transforms;  we 
have  defined 


Vi 


B 


<kl'j2)=^I  I  x(31'j2>  V”3 


•jlkl 


1  v° 


(15) 


which  requires  N  one-dimensional  transforms.  Thus,  N  +  N 

^  12 

one-dimensional  transforms  are  required  to  compute  the  single 

two-dimensional  transform. 

We  can  break  up  Equation  (13)  as  follows; 


A(ki#k2,k3) 


1 


n3-i 


I 


C(k  ,k  , j  )  w  ^3k3 

1  2  J3  3 


i3=° 


(16) 


which  requires  one-dimensional  transforms;  We  have  de¬ 

fined 


N  -1  N  -1 
1  2 


C(kl'k2'j3)  =  I  I  x(ji'32'i3)  w1"jlkl 


^2  :1=o  j2=o 


«2^2k2 


(17) 


which  requires  N  two-dimensional  transforms.  Thus,  N  N  one- 

^  1  z 

dimensional  transforms  and  two-dimensional  transforms  are 


needed  to  compute  the  single  three-dimensional  transform. 


/ 

/ 

/ 

/ 


; 

/ 

' 


4. 


ALGEBRAIC  DISCUSSION 


Equations  (9)  and  (10)  suggest  a  more  elegant  and  compact 
way  to  write  the  two  transforms.  We  define  the  vector  A  as 

the  transform  with  elements  (A)  =  A(k) ,  and  define  the  vector 

F  as  the  time  series  with  elements  (F)^  =  F(j)  .  The  process 

of  transforming  is  seen  to  be  equivalent  to  matrix  multipli¬ 
cation  by  a  matrix  W  whose  elements  are  (W)  =  w^k 


A  =  W+f 


(18) 


and  f  =  WA  , 


(19) 


where  the  dagger  indicates  Hermitian  conjugation.  Substituting 
(19)  into  (18)  gives  the  following  important  identity: 


+  + 

WW  =  W  W  =  I 


(20) 


This  is  the  definition  of  unitarity  for  the  transformation  W  . 
It  is  a  generalization  of  orthogonality  for  complex  matrices 
and  assures  Parseval 1 s  theorem: 


+  + 

A  A  =  f  f 


(21) 


W  preserves  "length"  between  the  two  domains.  The  identity  is 
actually  proved  by  writing  out  the  terms  in  the  product: 


N-l 


1  V  T  -|  3 m  r  -im^ 

N  L  l_exP(2TTi/N)  J  [exp(-2TTi/N)  J  =  6 


m=0 


7 


or 


N-l 


1 

N 


I 


m=0 


wm(j-k) 


(22) 


This  last  important  relation  is  seen  to  be  true  by  the  use  of 
a  phase  diagram: 


for  j  -  k  =  5  and  N  =  6  . 


^The  Cooley-Tukey  method  factors  the  W  matrix,  if  it  is 
a  power  of  two  in  order,  into  L  +  1  sparse  matrices,  where  L 
is  the  power  of  two: 


Multiplying  L  +  1  times  by  these  sparse  matrices  can  in  this 
case  reduce  the  computing  time  by  many  tens  of  times.  The 
factorization  is  proved  by  Good (4)  and  organized  for  computation 
by  Rader ( 3) . 


5.  HIGH-SPEED  CORRELATIONS  AND  CONVOLUTIONS 

* 

By  computing  Fourier  transforms  with  this  finite  Fourier 
series-like  method  an  important  condition  is  put  on  the  time 
series.  As  in  regular  Fourier  series  the  input  is  assumed  to 


-  8  - 


be  periodic  with  period  T  and  the  integrals  cr  sums  are  com¬ 
puted  over  a  single  period.  There  is  also  the  effect  of 
cutting  off  the  spectrum  at  the  folding  frequency.  Sines  and 
cosines  of  finite  wavelength  will  repeat  again  outside  the 
region  of  interest.  This  fact  in  itself  is  not  bothersome  but 
becomes  a  serious  complication  in  the  computation  of  convo¬ 
lutions  and  correlations.  Convolutions  and  correlations  as 
usually  computed  assume  the  time  series  to  be  zero  outside  the 
region  of  interest.  Therefore  the  integrals  or  sums  in  com¬ 
puting  them  are  summed  out  only  over  the  non- zero  region  of 
interest.  When  multiplying  together  two  finite  Fourier  trans¬ 
forms  (or  the  complex  conjugate  of  one  times  the  other)  the 
periodicity  of  the  time  series  means  that  elements  which  have 
been  shifted  past  the  end  of  a  period  reappear  at  the  be¬ 
ginning.  This  process  is  therefore  called  circular  convo¬ 
lution  or  correlation  and  its  effects  are  unavoidable  when 
straightforwardly  computing  lagged  products  with  finite  Fourier 
transforms.  This  is  illustrated  below: 


X1 

=  (3,  0, 

-1..  2) 

X2 

=  (-2,  2, 

-1,  2) 

(1. 

5,  3,  -1) 

for 

100% 

positive 

lags ; 

(1. 

-1,  3,  5) 

for 

100% 

negative 

lags . 

Circular  convolution  is  therefore  written: 

T-l 

Rc.(t)  =  Y  x,  (t)  x.(t  +  r)  (23) 

13  L\  1  j 

T"0 


9 


where  x^t  +  T)  =  xjt)  for  all  m  . 

The  proof  that  this  is  equal  to  the  transform  of  the 
product  of  the  two  finite  transforms  follows  below: 


T-l  T-l 

0,.  .  -tk 
ij  (t)  w  - 

I  l  Xi(T)  X,(t  +  T)  w'tk 

rt 

II 

O 

H 

II 

O 

l 

T-l  T-1+t 

I  I  xi 

(t)  x.(q)  W  (q“T)k  q  =  t  + 

li 

o 

►Q 

II 

T-l 

T-l 

7  X.(T)  Wtk 

L,  x 

T  =  0 

^  x.(q)  w  qk 

q=0 

=  A.  (k)  A  (k)  . 

/ 

I 

On  the  other  hand  the  transient  correlation  is  defined 
by  the  following: 


T-l-t 

Rij  (t)  =  I  Xi(T)  xj(t  +  t) 
T  =  0 


(25) 


Where  the  upper  limit  on  the  sum  simulates  the  desired  zeros 

in  the  time  series  outside  the  region  of  interest.  This  is 
lustrated  below: 


il- 
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Xj  =  (3, 

0, 

-1, 

2) 

X2  =  <-2 

,  2 

,  -1, 

3) 

*12  -  (1'  -4' 

6, 

-4) 

for 

100% 

positive 

lags; 

■  (1,  3, 

-3, 

9) 

for 

100% 

negative 

lags . 

n  c 

The  finite  Fourier  transform  of  this  R  is  thus  not  the  product 
of  the  two  individual  transforms.  However,  by  filling  zeros 
into  the  second  half  of  each  data  series  and  computing  their 
transforms  out  to  twice  their  actual  length,  a  good  estimate  of 
the  spectrum  may  be  obtained.  In  addition,  the  negative  lags 
in  the  correlation  appear,  thus  giving  a  more  mathematically 
satisfying  result.  This  is  illustrated  below: 

Xl  =  (3,  0,  -1,  2,  0,  0,  0,  0) 

X2  =  (-2,  2,  -1,  3,  0,  0,  0,  0) 

R-^2  =  ”4*  6,  -4,  0,  9,  -3,  3)  for  100%  positive  lags. 

The  two  modified  transforms  thus  are: 

2T-1 

F.(k)  =  V  X.(t)  w  tk  X.(t)  =  0,  T  <  t  <  2T-1 

1  4_i  1  1  —  — 

t=0 

2T-1  2T-1 

S.  . (k)  =  F.  (k)*  F.  (k)  =  V  X.(t)wtk  7  X.(t)  w"Tk 

i]  l  3  L>  i  L  j 

t=0  T=0 
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R'  (si 
ID 


*  )  F  .  ( k )  '  F 

L  i  j 


fk  )  wJ 


KS 


k=0 


2T-1  2T-1 


2m-l 


I  I  V:t)  VTl  I 


k ( t+S-T )  .  (26) 


w 


t=0  T=0 


k  =  0 


Now  from  (22)  the  last  sum  becomes  a  Kronecker  delta  function 
and  the  other  sum  is  collapsed  to  give; 


2T-1 


R  .  (s) 
iD 


-  I 


X  .  ( t )  X  .  ( t+s) 
i  D 


T 

R  .  (s 
iD 


t  =  0 


The  last  equality  following  from  the  original  assumption  that 

X . ( t )  =  0,  T  <  t  <  2T-1  .  Transient  correlations  for  100% 

i  _  _ 

lags  are  therefore  computed  by  forming  the  absolute  product  of 
two  transforms,  each  computed  out  to  twice  the  length  of  the 
original  data  series  with  zeros  filled  into  the  second  halves. 

Non-circular  or  transient  convolutions  are  computed  in 
much  the  same  way,  except  that  the  transforms  have  to  be  com¬ 
puted  out  to  a  length  equal  to  the  sum  of  the  lengths  of  the 
time  series  and  the  filter,  with  the  appropriate  number  of  zeros 
filled  into  each.  The  convolution  theorem  is  proved  in  the 
same  fashion.  . 

T+S- 1 

A(k)  =  a(r)  w  a(^)  =  0,  S<t<T+S-1 

t=0 
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X(t) 


0 


T<+<T+S-  1 


T+S-l 


X(k) 


•  I 


x(t)  w 


-tk 


t=0 


T+S-l  T+S-l  T+S-l  T+S-l 

A^k)  X(k)  wkU  =  £  £  a  (  t  )  X  ( t )  £  w  k(u-T_t) 

k=0  r«  0  t-0  k=0  j 

T+S-l  T+S-l 

J  A(k)  X(k)  Wku  =  J,  a (t )  xfu-t)  =  y(u)  (27) 

k  =0  t=o 


Where  y(u)  is  now  the  "filtered"  output  of  the  filter  a  acting 
on  X.  Convolutions  are  therefore  computed  by  forming  the  product 
of  the  two  transforms,  each  computed  out  to  a  length  equal  to 
their  sum  with  zeros  filled  into  the  extra  lengths.  Detailed 
procedures  for  these  computations  are  listed  in  Appendix  C. 
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APPENDIX  A  -  PROGRAM  LISTINGS 


FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


OOOOOOOOOOOOOOOOOOOOCSCJOOOOOOOO  nonooonnoooooooooooooooooo  ooo 


cv  it  f,e 

SUBROUTINE  UUUL(n()(,sIGN 


u 


HYPEH-RARID  FUURIFR  T  H ANSF  UHM  USING  CODlEY-TUKEY  ALGORITHM 


SEISMIC  DATA  LABORATORY#  ALEXANDRIA#  VA,  PROGRAMMED  ( 

26  FEBRUARY  1^66  BY  J»  F,  ClAfcRBOUT  tMIT)#  D,  R»  MCCOWAN, 

E.  A,  FLINN#  *ND  J.  GIBSON  (TELEDYNE) 

X  IS  A  COMPLtA  ARRAY  USED  FOR  THE  DATA  SERIES  AND  THE 

N 

TRANSFORM  •  1*E  NUMBEK  OF  ELEMENTS  OF  X  IS  L  *  2  « 

SIGN  •  -l.n  FOR  DIRECT  FOURIER  TRANSFORM  AND  *l,n  F OH  INVERSE 
FOURIER  TRANSFORM  <BUT  SEE  BELOW  FOR  ARRANGEMENT  OF  DATA  FOR 
INVERSE  IRANSFORM), 

FOH  DIRECT  TRANSFORM#  on  INPUI  THt  real  paht  of  x  contains  the 
data  SERIES  and  The  IMAGINARY  PART  of  X  is  ZERO,  ON  RETURN# 
the  Fourier  cusinf  series  expansiun  of  the  data  is  in  the  real 

PART  UF  X#  ANU  THE  FOURIER  SINE  SERIES  EXPANSION  IS  IN  THE 


N- 1 

IMAGINARY  pari  OF  X,  E ACM  CONTAINS  ONLY  2  ♦  1  NONHEDUNDANT 

POINTS,  The  COSINE  EXPANSION  is  SYMMETRIC  ABOUT  POINT  NUMBER 


2  1  ♦  1  And  Ire  SlNfc  TRANSFORM  is  antisymmetric  ABOUT 

THIS  POINT, 

FOR  EXAMPLE  •  N  *  3  ANU  D*T*  =  I  u  • » 1  •  *  I) •  *  0  •  *  0  • "  0  • »  0  •  *  0  •  *  • 

THEN  REAL  PARI  OF  X  «  t 0 . # 1 • # u . # 0 • • 0  • » 0 • » 0 • » 0 • 1  AND  IMAGINARY 
PART  OF  X  «  lu,#0#'0»'0*»Ut»0,»n*'0* >  INPUT, 

ON  RE  I  URN#  REAL  PART  UF  X  B  i 1 .  o  U  U#  •  7071»  0  •  *  *  •  '0  *  0  00* 

",7o7l#n.#«,.'u?l>  AND  IMAGINARY  PARI  OF  X  «  <o*»"*707l» 

■1  •  DO  U#  •  *  7  n7l#  n ,  #  ,  7  f)7l»  1*  u  D  U  #  •  7  0  7 1 1  •  POINT  NUMBER  l 
CORRESPONDS  IU  ZERO  FREOUENCY#  POINT  NUMBER  5  CORRESPONDS 
TO  PI#  THE  FULD1NG  FHEUUfcNCY,  ! 

TO  DO  AN  lNVtHSE  TRANSFORM,  T R£  CUSINE  AND  SINE  SERIES  MUST  BE 


N»  1 

FOLDED  over  about  POINT  number  i t  *1  BEFORE  CALLING 

COOL  RlTH  SIGN  «  ♦l,0«  SUBROUTINE  FTPACK  CAN  BE  USED  TO  DO 
THIS  FOR  YOU#  CONVERTING  AMPLITUDE  AND  PHASE  BACK  TO 
SINE  Aiyp  COSINg  IF  NbtD  BE. 


THERE  IS  A  SCALE  FACTOR  OF 


WHICH  COOL 


DOES  NOT  APPLY. 


THE  USER  CAN  APPLY  THE  SCALE  FACTOR  ElTHfcR  TO  THE  DIRECT  OR  TO 


THE  INVERSE  IRANSFORM# 
FOR  EXAMPLE#  GJVEn  THE 

call  cool<3»k#*i# 0> 

CALL  COOL<3*X#«l*o> 
would  chance  real  part 


-N/2 

OR  APPLY  A  SCALE  FACTOR  OF  2  To 
INPUT  DATA  AS  ABOVE#  THE  TWO  STATEMENTS 


of  x  to  ( q » #s# # o *  * o* * o *  *  o • *  • *  *  o » )  and 


c 

c 

c 

c 


c 

c 

c 
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IMAGINARY 


66 

PAH  I 


OF  X  TO  ^0*90t9099U**09*0,*Q,#Q,), 


dimension  xd)» iNTijei.acg) 
type  complex  x,u»w,hold 

EOUIVALENCF  <G,w> 

INITIALIZE 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 


IX  ■  2**N 

Pl2«6.?tJ<jlB5ij06 

FLX  «  LX 

FLXPI^«SIGNI*PIi!/FLX 
DO  lo  I  ■  1 »  N 
10  I  NT (  I )  «  2**<N-1  > 

LOOP  OVER  N  LAYERS 

DO  40  LAYEP  *  i»N 
NEfLOCK  ■  2»*ClAYER«i> 

LBLOCKeLX/NHLOCK 
LBHALF  *  LBLOCK/* 

START  SERIES  And  LOOP  OVfcK  BLOCKS  JN  EACH  LAYER 
NW  ■  o 

DO  A0  IBLOCK*!*  NBLOCk 
LSTART  ■  LRLOCK-T IBLOCK-i) 

COMPUTE  W  ■  LtXPl2.#PI*NW*SiGNi/LA) 

ARG»TLOAfF(NW)*»-LXPI? 

G  C 1  >  ■  COSFURGJ 

G<2>  «  SlNMARGJ 

ThjS  LAN  BE  SPEEDED  UP  bY  USJNq  A  TABLE  OF  COSINES 


COMI'UTe  ELEMtNTS  FOR  BOTH  HALF  S  OF  EACH  BLOCK 
DO  2u  I*X#LBHALF 

0  *  i+lstart 

K  ■  J+LBHALF 
Q  ■  X(K)*y 
X(K)  s  X< J)»Q 
X(J)  s  X(J)*U 
20  CONTINUE 


BUMP  UP  SERIES  BY  TWO  (NOT  ONE) 

00  32  I *2# N 

U  «  I 

LL *  1  NT  C I J ,ANU,NW 


THIS  LOGICAL  OPERATION  IS  A 
THE  APPROPRIATE  BJT  POSITION 
WORK  ON  IBM  FORTRAN  SYSTEMS, 


MASK  TO  DETECT  A  ONE  IN 
OF  NW,  THIS  STATEMENT  WILL  NOT 


IF  (LL)3i*3i«3q 


OOOOOOOOOOOO 


c: 

c 


u*  ?i 

.SUHHUUT  iwt  UtJCLCU,MMM,ltJ,,L,f  ,  ,  , 
lJ  I  MbNS  1  U'\i  rci»,xcz#i,.,  A.m*/> 

UIMt^lU.v  F(N»L>»X(^#ITfc!>1)»LAtt(l2/) 

MCLTICHANNfcL  CON  VOLM  I  J  ON  MOUMNfc  F  OH  I  APfcU  DAlM 

C:;  |'S  THfc  SUPSfcT  TAPfc  UF  U  A  f  A  CHANNfcLS 

‘  S"tt^T  UFE  OF  0  A  I  A  CHANELS 

f  is  .HHr  r  Titi  S"rp!xrirH  h,j,,*is  >ot<  ^ 

fr TAV  ct>NrMNlN&  A»  Lt AST  2.ITFJIT  POINTS 
HtSI  IS  l»,fc  nEXT  PnwtH  OF  ,wu  L AH|jtH  |HAN  LX*L 

L)  •  w  .  I’locnw  Aij  vjU^Y  1  VA6 


f<F:  W  1  im|.i  i  w  | 

Rpwimu  iu] 

*MU»  IFJI  JLAb 
N  S  L  A  rt  (  ?  ) 

L > =L  Ab ( J  i 
ISUM=LX»L 

l AH  I  3 i =u a- ( l-j  ) 

■^WITlrCTOI  )  |  am 

*-* 0  1  I  NU  =  1  ,  U 
I  I  PS  I  =  p  *  *  I  f.iu 
IF  *  I  SUM  -  I  If. Si 
?  NCClOL  a  J  NO 

do  lu  ;< 
l  Coni  iMit 
F’F'Iwi  iuuu.la.l 

luon  FOPMM,(&yHi«Atl  NtwS,  fcNH  Jk  I  .¥  LuoLCON, 
1X=  »Ifc,Un,  L=  ,1°) 

Slop 

3  Coni  i i\ i it 


UAlA  plus  F  I L  TfcP  Too  LONG  L 


I102MTfc.i1/2 
II  O^p^a  I  I  U?+<; 

IJO  lu  INS1,N 

call  tnAsfc(2*ntsi,x) 

'<t  AU  (  I  W  I  MvIi,HJ**is1,lvJ 

CO  11  TLS1.L 

11  ''I2»1L)sF(tn+(Il'*i)*m) 
call  cool (mcuoi.# a,-i.o > 
X(l*lJsX»l,i<«X(2,i )/|TfcS| 
X  (  2 > 1 )  s  u  *  ,(J 
IJO  20  IL  =  2,  1  |  ,12 


1“r-‘'^,‘*J-«I2.|TES,  ••  - .  TEST. 

20  X(l«iL)sSAvk 


-IL*2)**2-X<1, II  .)**2«X(2«TL>* 


30 

10 


Co  3U  IL=I  TU2H2*  ITFJjT 

X(l«iL)aA(^,i|psi»»lL  +  2) 
X(2» IL )=-X(2* I T  ts  T- 1 L+2  > 
CALL  COUL(MCOOL»X,*l0} 

WR 1 T t ( 1 0  I )(X<i»M>tM3L*LX) 
tML’  F  iLfe  IOT 
R  b  w I N  O  lur 
Rfcw  I  NO  I  .M  T 
Rfc  fUHN 


\ 


tNU 


nonnnnonnnnnnonnonnonnnnnnnn 
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SUBROUTINE  CUuLbH(N#X> 
DIMENSION  X ( 1) #  G ( 2) 
EQUI  VALt'^CF  <U*H» 


CUOLfcT-TUKfcY  ►0UH1ER  TRANSFORM  ON  REAL  T | ME  SERIES 
J.  F.  CLAbRBOUT  2«  JULY  1*66  ' 

INRUT  -  THE  heal  T I  Mb  SfcHIES  X ( % )#,.,# X ( L* ) 

OUTPUI  .  1 Hfc  UQMPLEX  FOURIER  TRANSFORM  X<l)#,.,,X(i  ♦  LX/2) 


NUTE  I  HAT  X  MUST  HE  I YPb  HEAL  IN  IHb  CALLING  PROGRAM#  AND 
DIMENSIONED  LX*i  THERb  (NUT  LX), 

SIZE  HER T H J C I  ION  -  LX  MUSI  Bb  L.b.  i6J«4 
(I.fc.,  N  MUST  bR  L.b.  1*) 


LX  I*H1*CK“1)*( J-i)/LX  ) 

X(J)  s  SUM  A‘K)*E 
K*1 

) ur  u  ■  l*LX'#;*i 
M 

ANU  WHCHF  LX  =  2 


ryPb  complfx  x, A»a,w,coNJG 

M  ■  N*>  j 

L  ■ 

CALL  COUL(M# X,-j.0  ) 

F  «  «i.i4159i;65»,)/fLOATF(LI 
W«X ( 1 ) 

X  ( 1 1  aRE AL  (  X  ( l )  )+A|MAa<X(i)  ) 

X ( L*1 ) sHbA|  (W)-AlMAG(H) 

LL  ■  L/2*l 
DO  lu  1 B2# LL 
J  «  L*I +4 

A«,5MCUNJR(X<  1  i  i*X(  J)  ) 

a«ts*(ru.')jR(x(  i  ( >«x(  j) ) 

U  «  I-i 
fl  ■  F*Z  1 
G(i)aCbSC ( F l ) 

G(2)«-SINF(Fi) 
d  *  d*u 

X(I)  a  A*(0,,i,|«B 
in  X(J)  ■  CUNjG( A)M Q,# 1, )*CQNJU(B> 
RETURN 
END 


OOOOOOOOOOOOOOOOOO 


OV  16  66 

subroutine  coolrlbr<n,x) 


THjs  COMPUTES  THE  HILBERT  TRANSFORM  OF  A  DATA  SERIES, 

THlISGRRSGRllMP^r!ASIfnF0UKIfcR  TRANSfwR«  ROUTINE  COOL 
TH|S  PROGRAM  IhANkS  TO  JON  CLAER80JT 


INPUTS  - 

N  *  tUG  <&ASb  2)  OF  NUMBER  OF  DATA  POINTS 

■  °AIA  seRIfS  TO  BE  TRANSFORMED 
I  M AG ( X  )  s  q 

OUTPUIS  * 

REAL(X)  a  x  AGAIN 
I  MAG  <  X )  a  HILBERT  TRANSFORM  OF  X 

THIS  calls  COUc 

dimension  xui 
TYPE  COMPLEX  X 

call  cool(n#X|*ji«q) 

M  *  2«*N 
»  M/**2 
00  1  I«Mi*m 
1  X(I)  *  <  u • # o • > 

X<1>  m  .B*X<1) 

X<Ml-i»  ■  ,5*XIMl#D 
CALL  COOL (N«X|*i«n) 

RETURN  n 

END 


J 


gv  it  tt 

subroutine  cooliwoin,x,sign,».wi 

c 

c 

C  This  USES  COUl  TO  COMPUTE  the  FOURIER  TRANSFORM  OF  TWO 

C  TIME  SERIES  AT  ONCE 

C 

C  INPUTS  • 

C  N  LOO  (BASE  2*  OF  NUMBER  OF  DATA  POINTS 

C  x  A  COMPLEX  ARRAY  OF  DATA,  The  first  time  series  Is  stored 

C  IN  THf  NEAL  PART  OF  X,  AND  THE  SECOND  IS  STORED  IN  THE 

C  IMAGINARY  PART  OF  X.  IN  OTHER  WORDS.  The  Two  SERIES  ARE 

C  MUUlPLExEU  |N  the  array  X, 

C  SIGN  •  *1.0  FOR  DIRECT  [RANSFORH,  THIS  SUBROUTINE 

c  HAS  NOT  been  CHECNEU  OUI  FOR  two  INVERSE  TRANSFORMS 

c  at  unce. 

c 

c 

C  OUTPUTS  • 

r  a  complex  Fourier  transform  of  the  first  data  series. 

C  I .1 ,,  THE  one  stored  in  the  real  part  of  X 

c  b  fourier  transform  of  the  second  data  series,  i.e..  the 

C  ONE  STORED  IN  THE  I HAb I  NARY  PARI  OF  X. 

C  BOTH  TRANSFONHS  ARE  OF  LENGTH  2«*lN-il  ♦  1  (SEE  COOL  WRITEUP) 

C 

c 

DIMENSION  X(ll<At,).6(il 

type  complfx  x.a.b.conju 

CALL  COOL(N.X|SIT*N> 

A  ( 1 »  R  ,»R(a(iI*CONJS(X(iM> 

6(j>  •  (u,,>.S)*(X(].)>CONJG(X(1)  )  ) 

Hb2**N 

00  \\i  K«*#H 

A(X)Rn,#R(XIK!‘GUNJS(X(H*2*Kl > 1 

10  B(XlR(0ig.«|).»Iw‘X(X>-CONJU(X(H*B.X)  It 

return 

END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


10 


20 


c 

c 

c 

c 

c 

c 


30 

35 

40 


C 

c 

c 


Subroutine  coqlvuivilx  x.if 

DIMENSION  Flli.XIJ.l, 

SINGLE -CHANNEL  CONVOLUTION  USING  COOL 


09  22  66 


iogetheh* ^and  transforhs^back **  ‘nd  FIl'ER'  multiplies 

INPUTS  . 

LX  LENUTh  OP  data 

PF  t f.NS™  0r  rILTER 

*  o{T.T/of°I;sF!0CiISTSxJ'J-„ONED  FILE.  IN  CALLING  PGM 

N  IS  T^u^ SMALLEST  number*  K,a  * 5  i  POWERH£p “*EE XCEBO I M« 

THE  SUBROUTINE  RETURNS  X  CONVOLVED  WlTu  r  „r 
if*LX’l.  STORED  CLOSE-PACKED  IN  x|!  '  °F  LENGTH 

23  SEPTEMBER  j?66  OHMCC 
CHECK  LENGTH  RESTRICTION 

nx*lf*lx 

°0  10  I«l#13 
N«2**J 

&Ti las 


error  Return 


lf«-lf 

return 


length  Of  Filtered  record  would  EXCEED  limit 


ncool  >! 

erase  working  space  in  X 
CALL  ERASE(N.NX,X<U41U 

multiplex  data  and  FILTER  in  X 
Co  30  1*1#LX 

J»LX-I4l 

X(1*J)  «  x ( J ) 

Co  3b  I»1,NX 

xcz,n  .  o.o 
Do  40  1*1* LF 
X(2,n  •  F ( I  ) 

TRANSFORM  and  FIDDLE 
fn«n 

CALL4COOLlNC0OL,X,.lt0I 
X,1#1>  «  X(l*i)*Xf2*i)/FN 


non  non 
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X(2,l)  •  0.0 
N2«N/2 

SSfca^MiassKsss . 

i  <4.*rm 

x!l.N2*l>*XCl»N2*i)*XC2»N2*i)/FN 

X(2»N2*1>*0*0 
N22»N2*2 
00  60  IL«N22#N 
XCl»!U**CliN»lL6|l 
60  X|2» Il>*  -XC2iN»lL*8l 

TRANSFORM  8AC*  \ 

CALL  COOLINCOOL»X|*1. 0) 

CL0Sfe*PACR  FILTERED  DATA-IN  X 


DO  70  I*l»  NX 
70  XII)  ■  X<1. D 
C 

e 

return 

END 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

CAUT 

C 

C 

C 

C 

C 

c 

c 

c 


U9  16  «6 

SUBROUTINE  FT2DUU0L<X»N»M,SlUNn 

two-dimensional  Fourier  thansform  using  cool 

INPUTS  . 

X(N,M>  ARHAy  To  BE  TRANSFORMED  IS  IN  REAL  PART  OF  X, 

*nd  the  imaginary  part  of  x  is  zero. 

N  FIRST  DIMENSION  OF  X 

M  SECOND  DIMENSION  OF  X 

SIGNI  ■  *1,0  ?!  OR  DIRECT  TRANSFORM#  FOR  INVERSE  TRANSFORM 

I0N  INVfcKSE  transform#  real  and  imaginary  parts  of  x  must 

BE  FOLDED  ABOUT  THEMIODlE  A$  JN  Cool.  SEE  COOL  WRITEUP, 

OOTPOI  - 

ON  RETURN#  RhAL  PART  OF  X  CONTAINS  COSINE  TRANSFORM.  IMAGINARY 
PART  OF  X  CONTAINS  SINE  TRANSFORM, 


D.  W,  mccowan 

DIMENSION  x(N#M! 
TYPE  COMPLEX  X 
FN  •>  N 
FM  •  M 


MAY#  1V66 


NCOOL  ■  L0GFTFNI/L0GFT2, T^l.E-6 
MCOOL  ■  L0RF(FM)/L0GF(?,)«i,E-6 
SN  *  l./SQRTF|FNI 
SM  ■  1,/SQpTF IFMI 
DO  l  IM*i,M 

CALL  COOL(NCOOL#X|i# im)#S1GN1> 

DO  1  IN  *  i#N 

1  X C IN# IM I  *  XC  JN# 1M)*SN 
CALL  MATHA63Cx»N.m,X) 

DO  2  1N«1*N 

CALL  COOL(  MCOOL#  X 1 1*  j  IK'*1)*M.  i>,  SIGNI ) 
DO  2  IM*1,M 

2  X(IN#1M)  ■  X< IN# IM)*SM 
CALL  MATRAA3(X#M#n.X) 

RETURN 

END 


oooonooooooonoo  non on 


«v  It 

SUBROUTINE  FTJDUUOKX.N.M.L.SIGNl  I 


THHEE-DIMENS1UNAL  FOURIER  TRANSFORM  USING  COOL 


INPUTS  . 

real  •'Art  o»  *  contains  the  three-dimensional  data  to 

TRANSFORMED*  The  IMAGINARY  PANT  OF  X  IS  ZERO, 


SE 


X  IS  DIMENSIONED  n  II  Ml  l 

SIGNS  «  FOR  DIRECT  TRANSFORM  AND  *!.„  FOR  INVERSE. 

TION  -  FOR  INVERSE  IRANSFORH*  DATA  MUST  BE  FOLDED  ABOUT  THE 

middle  as  jn  cool  see  cuol  rriteup. 

OU1PU I S  • 

ON  REfURN.  Real  PART  Of  X  CONTAINS  COSINE  TRANSFORM. 

AND  IMAGINARY  PART  OF  X  CONTAINS  SINE  TRANSFORM, 

D.  U.  MCCOVAN,  mAy.jVE* 

dimension  xir.m.li 
type  complex  x 

FL  ■  L 

LCOOL  ■  LOGFIFL.'LOGF is, l-i.E-6 
SL  •  1,/SQRTFIFLI 
DO  1  1L*1.L 

CALL  FY2UC00LIX(l,l, 1LI.N.M.BIGNI > 

1  Continue 

CALL  MATKA6G(X,N*M,L,X) 

DO  2  IN«1,n 
DO  2  |M"1,M 

CALL  COOLILCOOL.  A{  j*(  IN-iI*L*HM'i  I-LaN.  i  • ,  >  , sigri  I 
DO  2  IL*1.L 

2  X(iN.IM.iL)  •  XUN.1M,  IL1-SL 
CALL  MATNA6S(X,L.N*M,XI 

return 

End 


o  o  o 


subroutine  ma^ra© j (  a, n, m,  b  > 
dimension  a < 2| i > » b ( 2< i ) 

MATRIX  TRAnSHUSE  ON  COMPLEX  ARRATS 

maski*oouooooooouooooib 

MASK2.77/777777//77776B 

nm°n*m 

DO  1U  1*1, NM 
B(l»  I )  ■  A  ( 1,  I ),0R»MASK1 
10  B(2, 1 )«A<2,  I I 
JF*0 

ASSIGN  3U  TO  KShH 
DO  100  1*1, NM 
GO  TO  KSRH,  ( 3q,  S>U  } 

30  JF  =  JFu 

UL=B(l, JF ) .AND.MASKl 
IF  <LL)3Q*3o»«o 
40  JO* JF  «l 

ASSIGN  50  TO  KSWH 
TEMPHl.b(i.jF) 

TEMPB2*B(2, JF) 

50  Jl»JO/N*XMOUF( JO*N)*M*i 
T£MPAl»B  C 1, J1 ) 

TEMPA2*b ( 2, J1 ) 

BCl,  Jl)*rEMPBi,AND,MASK2 

9(2, Jl)* TEMPB2 

TEMPBl«TtMpAl 

TEMPB2.TEMPA2 

J0*J1*1 

IF ( Jl-JF ) 6  o , 6o,lUO 
60  ASSIGN  30  TO  KSWH 
100  CONTINUE 
RETURN 
END 


onooooooonooooo 


c 

c 

c 

c 

c 
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subroutine  specimumut.  jt.mt.x.l.lf.s» 
dimension  x<2(i>»s<2,i> 

DIMENSION  x(2fN«N|Uf||X( 2*2* LX* 2 1  •  SC  2#LF ) 


SPECTRAL  M1TMIX  FOR  T*PED  D*T* 


SS:i.7U8T.2.t  MME*  0F  TN°  |N  LfeNfiTH  AN0  0N  TApE  *T  Ik  SUB  RET 
FORMA! •  SREClR&L  MATHJX  IS  RbTUKNED  As  A  F6J  COMPLEX  MATRIX 
1 N  X  i 


it-inmut  subset  tape 
jt -scratch 
kt-scnatcm 

X-MONKInG  ARNAy  And  RETURNED  SPECTRAL  MATRIX 
L'NUMMER  of  TIMES  TO  SmOOIm 
LA •RETURNED  LfcNGTW  OF  SPECTRAL  ESTIMATES 
S*hOHK|nG  ARM»t 


PHOSNAm  TOO  COMPLICATED  TO  DESCRIBE,. , 

REMIND  IT 
remind  or 
remind  mi 
readiitilost.n.u* 

LX2»2*L* 

NCODLnLUUF (f LOA I n LX  I  T JLUfiF I m. n »♦, . „E.6 

N$0«N»N  " 

IF*LX/2»»L*1 

lX2P2aLXB+p 

LXA»A*LX 

LX2P2T2«BalX2P2 

LXPl»LX»l 

LX2PS>LX2«3 

IDC«o 

LF2»2»LK 

»RITE(JTIL0S1,N.LX2P2 
HR1TECKTILOS1|N#LX2Pp 
DO  10  IN*l,N 
CALL  ERASElLXA.XI 
READUTMXIl.MT.HfL.LXI 
CALL  COUUNCUOLAl.X#. l.iil 
»RITEIjTMX(M|.M*1,LXjP2l 
MRITE(MTHxlMI,M«,,LX2P2l 
1(  CONTINUt 
end  file  jt 

END  FILE  KT 
REMIND  IT 

remind  ji 

REMIND  Ml 
DO  1  IN*1.N 
INDe1N»1 

CALL  SKIMRECIIN,MT) 
«E»OI«TFlX(MI,M»l,LX2PF| 

call  DOTtM(x»x»L*pi,i((LX*pji> 

CALL  SM0UTmIX(LX2P3),lxPi.LI 
CALL  DISC63I IDC» 1,X(LX2PJ) ,LF2I 
IDC»IDC«» 

DO  6y  JN>1ND.N 
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«t*D<KT )  (X(M)#M*LX2P3#LXk:Pi;Ti:) 

^ *LJ“  ^It”(X*K“-X2P3>»LXP1,)((LX2P3?) 
CALL  SHOUTwCXJLX^pjJ.lxP!^) 

Call  0ISC63< IDC#i,x<LX2P.J>#L>2> 

IDC* IDC*V  * 

6q  CONTINUE 
HfcWJNO  ftl 
IS*Vfc»KT 
KT«JT 
JT«lSAVb 
1  CONTlNUb 
I  DC*  y 

DO  2t>  I Nrl» N 
1 ND« 1 N*i 

CALL  CISCfcSCIDL^g^iLF?) 

IDC« 1 DC*V 


INDEX«IN*< IN-i)«N 
DO  2«>  IL*1,LF 
*<1» lNDbX)*S(i# 1L| 

X(2#  I Np&X > «s < 2#  1L» 
INOEXs]'NUEx*NSQ 
26  CONTlNUb 

DO  2'/  JN*1nO«N 

CALL  DISC63(  mc#mS#LF2> 
1 D  C  *  I CC* V  L  * 


2B 

n 

25 


lNDEXi«IN*(UN*i)*|g 

INDEX2bUN«v IN»? I*N 
DO  2«  IL=i,LF 
X(i»INDEXi>«S(i#1L) 

X  C  2* lNDbXi)*S(2# 1 L  > 
X(l#  INDEXgjeSCj,#  lL> 
X(2# 1NDEX2)««S(2»  J L I 
INDEXjbInDEXi^NUU 
INDEX2«JNDFX2*NSU 
CONTINUE 

continue 

continue 

rituhn 

END 


OUU 


QV  16  6<> 

SUBRUUT  1  Nfc  UUTEMtx,  Y,UZ> 

DIMENSION  XC2*D*V<2»D#Z(2#U 
DO  1  R«i,l 

SAVENMCl# IL>*Ytl* !L>*X<2» 1D*Y<2, ID 
SAVEI*X(l#tU*YU,  Il>-V<2,  R)*YC1,  R) 
Z ( 1  *  I L ) »SA vfcH 
1  Z(2»  IL)»6Avbl 
HE  TOWN 
EM) 


SUBROUTINE  SM0OTHjx,LENOTHR) 

IhIS  HANNlNU  HQUTJNE  TMANNS  TO  J  C'LAfcHUOUT 

DIMENSION  X ( 2i LbNQTH ) 

U*LENGTH 
LFMjsLF-1 
DO  j  IL*i»L 

*<  1*  1)«U*5*X<  1>*0,‘>*X(1»2) 

X  (  2» 1 ) « u«  o 

X(  1»LF  )«0.5*X  I  l#Ur  >*f,,5*X(  i»LF-i  > 

X( 2»LF  >*0t  n  1 

I  NO*  2 

DO  2  Jl^RFM*,* 

X( 2# Ju>«u. 26*x«  2* JL-i)*u.5*X»  2# JL)*0, 2&*X( 2# JU+1 ) 

X(2» 1ND»»X(2»JL) 

IND  =  1ND«>1 
2  CONTlNUfc 

X(i»  JN0)«X<iRF  I 
X<2»  1ND»«X(2RF  ) 

LF'LF/2*1 

IFMjRF-l  i 

i  continue 

RE  I  URN 

end 


Ooooooooooooo 


subroutine  disca^i jblock* iswitcu.v  km 

DIMENSION  y (N)  -,1d,*0C''»‘swiTCH#X,N) 

ii "«?Ki!"eK.uTis  tT;&,n  ~ 

»« in 

SW  TCH.O  GIVES  A  HEAD  FROM  THE  DISC 

x  IS  Tw'JSEMjJiMS*  "R1TE  °N  THE 
n  is  the  number  of  words  to  transfer 
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c  ,H‘S  ,<0Ur‘NE  "USI  6E  SUPP‘-*t‘’  SX  'Ht  USER  OR  INCLUDED  siRRRY 

return*  ********************************************************* 


o  o  o 
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SUBROUTINE  ERASbtN.X) 
DIMENSION  X(Nt 

ERASE  N  WORDS  IN  X 

DO  1  I*1#N 
1  X  C 1 >  *0 . U 
RETURN 
END 


o  o  o 


SUBHOUTlGE  SK]PHtC(N#|TAB£i 

SMP  H  LOGICAL  RECOROS  „„  T,pt  lJApg 

00  1  I«1»N 
1  «EAD(ITAHE»L0ST 
hetunn 
end 


APPENDIX  B  -  PROGRAM  WRITE-UPS 


FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 


DIGITAL  COMPUTING  SECTION 
A.  IDENTIFICATION 

T_it_le:  Hyper-Rapid  Specialized  Cooley-Tukey  Fourier 
Transform 

COOP  Identification;  G612-COOL 
Category;  Fourier  Transform 

Programers;  J.  F.  Claerbout,  D.  W.  McCowan ,  J.  L.  Gibson, 
and  E.  A.  Flinn 

Date;  26  February  1966 
B „  PURPOSE 

To  compute  the  Fourier  series  expansion  of  a  real-or 
complex-valued  date  series,  or  the  data  series  from  the 
complex-valued  Fourier  series  expansion. 

C .  USAGE 

1 •  Operational  Procedure  and  Parameters; 

This  is  a  CODAP  subroutine  with  a  FORTRAN-63  calling 
sequence  CALL  COOL  (N,  X,  SIGN).  X  is  a  complex  array 
used  for  the  data  series  and  the  transform;  the  number 
of  elements  of  X  is  L  =  2N?  SIGN  *  -1.0  for  a  direct 
Fourier  transform,  and  +1.0  for  an  inverse  Fourier 
transform  (but  see  below  for  arrangement  of  data) . 

For  the  direct  transform;  on  input  the  real  part  of 

X  contains  the  data  series  and  the  imaginary  part  of  X 

is  zero.  On  return,  the  Fourier  cosine  series  expansion 

is  in  the  real  part  of  X,  and  the  Fourier  sine  series 

expansion  is  in  the  imaginary  part  of  X.  Each  contains 
N-l 

only  2  +1  non  redundant  points;  the  cosine  expansion 

.  N-l 

is  symmetric  about  point  number  2  +1  and  the  sine 
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transform  is  antisymmetric  about  this  point. 


For  example:  N  =  3  and  data  =  (0.,  1.,  0.,  0., 

0.,  0.,  0.,  0.);  Re (X)  =  (0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,) 
Im(x)  =  (0,,  0.,  0.,  0.,  0.,  0.,  0.,  0.,)  on  input.  On 
return,  Re(X)  =  (1.000,  .7071,  0.,  -.7071,  -1.000, 

-.7071,  0.,  -.7071);  Im(X)  =  (o.,  -.7071,  -1.000, 

-.7071,  0.,  .7071,  1.000,  .7071).  Point  number  1  corre¬ 

sponds  to  zero  frequency;  point  number  5  corresponds  to  n 

For  inverse  transform:  the  cosine  and  sine  series 

N-l 

must  be  folded  over  about  point  number  2  +1  before 

calling  COOL  with  SIGN  =  +1.0 

-N 

There  is  a  scale  factor  of  2  which  COOL  does  not 
apply.  The  user  can  choose  to  apply  the  scale  factor 

•  s 

either  to  the  direct  or  to  the  inverse  transform,  or  to 

-N/2 

apply  a  factor  of  2  to  both.  For  example,  if  COOL 

were  called  with  the  transform  example  above,  the  result 
would  be  Re(X)  =  (0.,  8.,  0.,  0.,  0.,  0.,  0.)  and 
Im (X )  =  (0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.)  . 

3..  Space  Required:  Approximately  200^Q  exclusive  of  X. 

The  largest  series  that  can  be  transformed  in  a  32K 
core  machine  is  8K. 

4.  Temporary  Storage  Required:  None.  Other  versions  of 
this  program  have  an  auxiliary  storage  for  the  cosine 
table  and/or  a  table  of  bit-reversed  numbers.  COOL 
computes  its  sines  and  cosines  as  it  goes,  and  uses  an 
algorithm  due  to  J.  F.  Claerbout  for  calculating  the 
bit-reversed  numbers. 

5.  Printout :  None. 

6\  Error  Printouts:  None. 
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7,  Error  Stops;  None. 


8*  Input  and  Output  Tape  Mountings-;  Net  Applicable, 

9*  Input  and  Output  Forma, s;  Not  Applicable. 

10.  Selective  Jumps  and  Stops.  None. 

11-  Timing;  Time  is  proportional  to  N’2N  .  Transforming 
8192  on  the  CDC  1604-B  requires  25.0  seconds. 

12  Accuracy:  Calling  COOL  returns  the  original  to  about 
nine  decimal  places. 

1  '  ^autions  to  User:  See  Operational  Procedure  above. 

14-  Configuration:  Standard  COOP. 

*  References  J.  W.  Cooley,  1964  "Harm  -  Harmonic  Analy¬ 
sis;  Calculation  of  Complex  Fourier  Series*:  IBM 
Watson  Research  Center,  Yorktown  Height,  New  York. 

J.  W.  Cooley  and  J.  W.  Tukey,  1965,  An 
Algorithm  for  the  Machine  Calculation  of  Complex 
Fourier  Seriess  Math,  of  Comp.,  Vol .  19,  pp ,  297-301. 

Cc  M.  Rader,  1965,  "Algorithm  for  Rapid 
Digital  Computation  of  a  Spectrum,  MTT  Lincoln  Labo¬ 
ratory,  Lexington,  Massachusetts. 

D.  W.  McCowan.  1966,,  "Practical  and  Compu¬ 
tational  Aspects  of  Fourier  Transforms,"  Teledyne,  Inc  , 
Alexandria,  Virginia. 


Writeups  of  the  following  SDL  programs 


C00LTW0: 
COOLEST  s 

COOL EXT : 

FT  3DC00L : 
FTPACK: 


Does  two  Fourier  transforms  at  once. 
Does  Fourier  transform  of  data  series 
of  length  other  than  a  power  of  2„ 
Does  Fourier  transform  of  16384  data 
points . 

Three-dimensional  Fourier  transform 
Driver  for  COOL  (converts  amplitude  and 
phase  to  sine  and  cosine,  does  the 
folding,  etc.) 
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D .  METHOD 


N 

Given  a  time  series  X(l),  1,  L  (where  L  *  2  )  assumed 
to  be  periodic  outside  the  given  range,  COOL  constructs 

N-l 

Y(K)  “  SUM  X(J)*WJK  K  =  0,,  L  -  1 

J*0 


where  W  =  exp  (-2  i/L)  for  time-frequency  transform,  and 

W  =  exp  (+2  i/L)  for  f requency-t ime  transform.  The  algo- 

N 

rithm  is  efficient,  requiring  N* 2  multiplications  rather 

,  _  2N 

than  2 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA ,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A .  IDENTIFICATION 

Two  and  Three  Dimensional  Fourier  Transform  tdckage 
COOP  Identification:  G615  FT2DC00L,  FT3DC00L 
Category:  G6  Time  Senes  Analysis 

Programer:  D.  W.  McCowan 

Date:  20  April  1966 

B .  PURPOSE 

The  subroutines  in  this  package  compute  two  and  three 
dimensional  Fourier  transforms.  Their  names  are:  FT2DCOOL , 
FT3DCOOL,  COOL,  MATRA63,  and  SCALE.  As  with  COOL,  the 
dimensions  on  the  data  must  be  a  power  of  two. 

C o  USAGE 

1.  Calling  Sequence: 

CALL  FT2DCOOL  (X,  N,  M  SIGNI) 

and 

CALL  FT 3DCOOL  IX,  N,  M,  L,  SIGNI) 

2 .  Arguments : 

X,  the  complex  array  in  which  the  data  is  supplied  and 
in  which  the  Fourier  transform  is  returned.  If 
real  data  is  supplied,  it  must  be  put  into  the  real 
part  of  X  and  the  imaginary  part  must  be  erased. 

N,M,L,  the  dimensions  of  X.  Each  of  these  numbers  must  be 
a  power  of  two.  The  number  of  complex  points  in  the 
Fourier  transform  will  be  N/2  +  1,  M/2  +  1,  and 
L/2  +  1  in  each  direction. 

a  switch  determining  the  type  of  transform  to  be 
performed.  SIGNI  =  -1.0  gives  an  direct  transform 
(time  to  frequency),  and  SIGNI  =  +1.0  gives  the 
inverse . 
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SIGNI , 


3. 


Space  Required j  500  locations 


4.  Temporary  Storage?  None 

5.  Alarms  and  PrintouLs  None 

6.  Error  Returns  None 

7.  Error  Stops?  None 

8.  Tape  Mountings:  None 

9.  Format  s :  None 

10.  Jumps  and  Stop  Settings:  None 

11  •  Time  Required:  Three-d?.mensional  Fourier  transforms 

require  NM  +  NL  +  ML  one-dimensional  Fourier  transforms „ 
Two-dimensional  Fourier  transforms  require  N  +  M  one¬ 
dimensional  Fourier  transforms.  For  the  timing  of  one¬ 
dimensional  Fourier  transforms,  see  References. 

12 .  Accuracy:  Same  as  COOL. 

13.  Cautions  to  Users:  None 

14.  Equipment  Configuration  Standard  COOP 

15 .  References.  Writeup  of  UES  G612  COOL  3/30/66 
D .  METHOD 

The  direct  2  and  3-dimensional  Fourier  transforms  are 
defined  as: 


A  (  j  i  '  j  2  ) 


N-l 


and 


M-l 

X(k  k  )  W  -jlkl  W  ~j2k2 
L  Z  L  z 


N-l  M-l  L-l 


A(il'j2’j3)  *  III  ^I'W 


NML 


k  =0  k  =0  k  =0 
12  3 


W1*Jlkl  »2-32ki  w  -J3k3 
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Where  =  exp(2^/N);  *2  =  expUm/N);  W3  =  exp(2m/L)  . 

The  two-dimensional  transform  is  broken  up  onto  N  +  M 
one -dimensional  transforms  and  the  three-dimensional  trans¬ 
form  is  broken  up  into  L  two-dimensional  transforms  and  NM 
one-dimensional  transforms. 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A.  IDENTIFICATION 

Title:  Fourier  Transform  of  Two  Data  Series  Simultaneously 

COOP  Identification:  COOLTWO 

Category:  G6  Time  Series  Analysis 

Programer:  E.  A.  Flinn 

Date;  10  June  1966 

B .  PURPOSE 


To  compute  the  Fourier  series  expansion,  using  COOL 
(q.v.),  of  two  data  series  simultaneously. 

C .  USAGE 

t 

1*  Operational  Procedure:  This  is  a  FORTRAN-63  subroutine 

with  calling  sequence 

CALL  COOLTWO  (N,  X,  SIGN,  A,  B). 

2 .  Parameters : 

N  is  the  log  (base  2)  of  the  number  of  elements  in  X; 

X  contains  the  two  data  series,  multiplexed  in  one 
complex  array,  so  that  Re(X)  contains  one  series  and 
Im(X)  contains  the  other; 

SIGN  —  -1.0  .  The  program  has  not  yet  been  checked 
out  for  inverse  transformation; 

A  is  the  complex  (cosine  and  sine)  transform  of  the  data 
series  stored  in  the  real  part  of  X; 

B  is  the  complex  Fourier  transform  of  the  data  series 
stored  in  the  imaginary  part  of  X; 

A  and  B  are  both  of  length  2** (N  -  1)  +  1  , 

3-  Space  Required:  about  70, ^  excluding  arrays. 
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^ ■  Temporary  Storage  Requirements;  None 

5.  Printouts ;  None 

6*  Error  Printouts;  None 

1 •  Error  Stops;  None 

8.  Input  and  Output  Tape  Mountings:  None 

9.  Input  and  Output  Formats:  Not  Applicable 

10,  Selective  Jump  and  Stop  Settings:  Not  Applicable 

n*  Timing  is  proportional  to  N‘2N;  transforming 

8192  data  points  on  the  CDC  1604-B  requires  25.0  seconds. 

12.  Accuracy:  Same  as  COOL 

13 *  — a-U— — ns  to  User;  This  program  has  not  been  checked  out 

for  inverse  transformation.  This  program  does  not 

apply  the  scale  factor  2_N,  since  some  users  may  wish 
to  apply  the  scale  factor  to  the  inverse,  rather  than 

the  direct  transform.  The  number  of  data  Points  must 
be  a  power  of  2 . 

14.  Configuration:  Standard  COOP 

15.  References:  Writeup  of  UES  G612  COOL 
D„  METHOD 

The  method  is  due  to  J .  W.  Cooley  (see  Reference  2  in 
main  body  of  this  report) 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 


A.  IDENTIFICATION 

Title ;  Spectral  Matrix  Estimates 
COOP  Identification;  G618  SPECTRUM 
Category;  Time  Series  Analysis 
Programer:  D.  W.  McCowan 

Date;  10  July  1966 

B .  PURPOSE 


This  is  a  package  of  three  FORTRAN-63  subroutines  for 
computing  an  estimate  of  the  spectral  matrix  for  N  channels 
of  data  stored  on  magnetic  tape,.  It  uses  the  hyper-rapid 
Fourier  transform  routine  COOL,  and  makes  use  of  two  tapes 
and  the  disc  to  cut  running  time  to  a  minimum,.  The  names 
of  the  three  routines  in  the  package  ares  SPECTRUM,  DOTEM, 
and  SMOOTH.  In  addition  to  these,  three  more  subroutines 
are  assumed  to  be  on  the  system  tape;  they  are;  COOL, 
SKIPREC,  and  ERASE.  Since  all  other  routines  are  called 
internally  by  SPECTRUM,  only  the  calling  sequence  for  it 
will  be  given. 

C .  USAGE 

1.  Calling  Sequence; 

Call  SPECTRUM  (IT,  JT,  KT ,  S,  NS,  LF ,  X) 

2 .  Arguments; 

IT,  the  input  subset  tape  number  on  which  the  N  channels 
of  data  are  written.  The  length  of  eaqh  channel 
must  be  'exactly  a  power  of  two. 

JT,  the  number  of  a  scratch  tape. 

KT,  the  number  of  a  scratch  tape. 
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a  triply  subscripted  FORTRAN-63  complex  array  used  both 
for  internal  manipulation  and  to  return  the  computed 
spectral  matrix  as  a  N  by  N  by  LF  complex  array  with 
subscripts  varying  in  that  order.  Here  N  is  the  number 
of  channels  read  from  the  input  tape  label  and  LF  is  the 
smoothed  length  of  each  spectral  estimate.  This  array 
must  also  be  4*LX+4  locations  in  length,  since  it  is 
also  used  for  internal  computations,,  LX  is  the  length 
of  the  input  data  channels  read  from  the  input  tape 
label.  Remembering  that  there  are  two  locations  used 
for  each  complex  number,  the  total  dimensions  on  S  in 
the  main  program  must  be  2*N*N*LF  or  4*LX+4,  whichever 
is  the  larger.  It  is  usually  convenient  to  dimension 
it  as  a  complex  N  by  N  by  L  F63  array  in  order  to  fa¬ 
cilitate  use.  L  here  is  a  number  chosen  so  that  S  will 
be  large  enough  as  described  above. 

the  number  of  times  to  apply  the  hanning  smoothing 
operation  to  the  original  estimates. 

the  returned  length  of  the  spectral  estimates.  This  is 
computed  from  the  formula: 

LF  =  (LX/(2**NS)  +  1 
LF  must  not  be  larger  than  129  . 


an  array  used  for  internal  manipulation,  containing  at 
least  2'*LF  locations. 

Space  Required:  502  locations 

Temporary  Locations;  None 

Alarms  or  Special  Printout:  None 

Error  Returns:  None 

Error  Stops:  The  subroutines  stops  if  length  of  filter 
plus  length  of  and  chcinnel  exceeds  213 

Tape  Mountings.  See  Arguments 

Input  and  Output  Formacs:  See  Arguments 


Jump  Settings s  None 


11.  Time  Required.  A  10-channel,  4096-point  NS  =  6  case 

takes  approximately  10  minutes  of  1604 
time . 

12.  Accuracy:  Single  precision 

13.  Caution  to  tsars  The  subroutine  as  written  requires 

that  the  data  series  should  contain 
a  number  of  points  exactly  a  power 
of  two . 

14.  Equipment  Configuration:  Standard  COOP 

15.  References ;  Writeup  of  subroutine  UES  G612  COOL  6/1/66 

Writeup  of  program  UES  Z24  SUBSET 

Stockham.  T.  G.,  1966,  High  Speed  Convolution 
and  Correlation,  AFIPS  Proceedings 


METHOD 


The  spectral  matrix  elements  S. , (k)  are  usually  defined 

ID 

as  Fourier  transforms  of  correlation  functions  R  .(t).  How- 

-  D 

ever,  it  must  be  realized  that  these  correlations  are  transient 
correlations  where  the  functions  are  considered  to  be  zero 
outside  the  reqicn  of  interest  and  100%  lags  are  taken.  They 
are  defined  as  follows 


T-l-t 

R.  t,  »  V  x  i t  x.(r  +  t' 
1]  -  i  3 

t=0 


T-l 


R  =  V  * 

1  j  U  1 


x  .  (r  -  t  j  =  R.  . 

3  31 


(t; 


T*t 


The  spectral  matrix  element  is  then 


T-l  T-l 


T-l  T-l-t 


3  (kl  -  Y  V  x. (t)  x . it  -  t  W  +  Y  Y  x  ( rl  x .  ( r+t ) 
13  -  1  3  *  1,  u  1  3  . 


t=0  —t 


t=l  T=0 


-tk 


Si, 
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This  can  be  shown  to  be  equivalent  to; 


sij(k)  =  Vk)  Fj(k)' 

* 

where 


T~1 

f,  (k)  =  y  x,  (t)  w'~ 

i  L,  i  2 

t=0 


This  is  recognized  as  the  Fourier  transform  of  the  input 
data  computed  over  twice  its  length  with  zeros  filled  into 
the  second  half.  The  Cooley-Tukey  hyper-rapid  Fourier 
transform  routine  COOL  is  used  to  provide  the  high  speed 
necessary  here. 

Each  spectral  matrix  element  is  originally  T  +  1 
complex  points  long  between  DC  and  the  folding  frequency. 
It  is  then  smoothed  with  a  banning  window  NS  times  to  its 
final  length  of  LF  points. 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A-  IDENTIFICATION 

Title:  Hyper-Rapid  Specialized  Cooley-Tukey  Fourier  Transform 

(direct  only) 

COOP  Identification:  G617-COOLER 

Category:  Fourier  Transform 

Programer:  J.  F.  Claerbout 

Date:  27  July  1966 

B.  PURPOSE 

To  compute  the  Fourier  series  expansion  of  a  real-valued 

time  series. 

C  .  USAGE 

1.  Operational  Procedure:  This  is  a  FORTRAN-63  subroutine, 
with  calling  sequence  CALL  COOLER(N,X).  This  subroutine 
calls  COOL . 

2.  Parameters ?  On  input,  X  is  a  real-valued  time  series 
containing  LX  points,  where  LX  =  2^,  N  is  restricted 
to  be  14  or  less.  On  return,  X  contains  ^LX+1  complex 
points  of  the  Fourier  transform  of  the  data,  with  the 
real  and  imaginary  parts  multiplexed  together  -  i.e., 
on  return  X  can  be  thought  of  as  a  complex  array,  with 
the  cosine  transform  in  the  real  part  and  the  sine 
transform  in  the  imaginary  part. 

X  must  be  dimensioned  at  least  LX+2  in  the 
calling  program,  (i.e.,  ^LX+l  complex  points' 

3.  Space  Required:  very  little  , 

)  ‘ 

4.  Temporary  Storage  Required:  none 
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5, 


Printout:  none 


6 .  Error  Printouts;  none 

7.  Error  Stops:  none 

8.  Input  and  Output  Tape  Mountings:  not  applicable 

9.  Input  and  Output  Formats:  none 

10.  Selective  Jumps  and  Stops:  none 

N 

11.  Timing:  Time  is  proportional  to  N .2  ;  transforming 
16384  points  on  the  CDC  1604-B  requires  45.0  seconds. 

12.  Accuracy:  About  nine  decimal  places 

13.  Cautions  to  User:  On  return,  the  real  and  imaginary 
parts  of  the  transform  are  multiplexed  together.  X  must 
be  dimensioned  at  least  LX+2  in  the  calling  program,  not 
LX.  This  subroutine  will  not  do  an  inverse  transform. 

14.  References :  Writeup  of  UES  G612  COOL. 
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SEISMIC  DATA  LABORATORY 
AL EXANDRIA ,  VT  RG1N1A 

DIGITAL  COMPUTING  SECTION 


A.  IDENTIFICATION 

Title:  Multichannel  convolution  in  the  frequency  domain, 

for  taped  data. 

COOP  Identification  L'fiS  G620  COOLCON 
Category.  G6  Time  Series  Analysis 
Prog  r  a  me  r :  D.  W.  McCowan 
Da te •  22  September  1966 

B .  PURPOSE 

This  subroutine  convolves  data  channels  on  the  input 
subset  tape  with  a  multichannel  filter  stored  in  core, 
working  entirely  in  the  frequency  domain.  The  result  is 
written  in  subset  format  on  another  tape. 

C .  USAGE 

1.  Operational  Procedure  This  is  a  FCRTRAN-63  subroutine 
with  calling  sequence. 

CALI.  COOLCON  (I  NT,  IOTL.F,Xl 

2 .  Parameters 

INT  is  the  number  of  the  inpuc  tape  unit . 

IOT  is  the  number  of  the  output  tape  unit. 

L  is  the  number  of  points  in  the  filter  (see  restriction 

below) . 

F  is  the  multichannel  filter,  dimensioned  F{N,L)  in  the 
calling  program,  where  N  is  the  number  of  channels  on 
the  input  subset  tape. 

X  is  a  working  array,  dimensioned  X(2,IT)  in  the  calling 
program,  where  IT  is  the  least  power  of  2  such  that 

IT 

2  >  L  +  LX 

where  LX  is  the  number  of  data  points  in  the  input 
channels . 
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Restrict  ion  on  length  of  data  and  length  of  filter. 

LX  +  L  must  not  be  greater  than  213  (8K> . 

3.  Space_Required:  Very  little  in  addition  to  arrays. 

4'  — P-rary  Stora96  Required  2‘2IT  working  space  plus 

12710  for  the  subset  tape  label. 

5-  Printout ;  None. 

6.  Error  Printouts;  rf  r.+rv^i^l  t _ 

- =.  it  these  numbers  are  printed 

w’-th  an  error  message. 

7.  Error_  Stops;  If  L+LX>213,  the  subroutine  stops  the 

calling  pr-gram. 

8'  — ^  wd  Output  Tape  Mountings;  See  Parameters  above. 

jnput  and  Output  Formats.  Compatible  with  UES  Subset 
(See  Writeup) . 

10-  Selective  Jump  and  Stop  Settings;  None. 

n.  Timing;  Dominated  by  two  Fourier  transforms  using  COOL 
f°r  each  channel  to  be  filtered.  The  length  of  trans¬ 
form  is  2  (See  Writeup  of  COOL). 

Accura .y.  This  yields  the  same  numbers,  to  ten  decimal 
places  which  would  be  computed  by  convolving  the 
filter  and  data  series  in  the  usual  way. 

13.  Cautions  to  User;  None. 

34-  Configuration :  Standard  COOP. 

15.  References:  Writeups  of  UES  0612  COOL,  UES  Z24  SUBSET, 
and  UES  G617  COOLER. 

D,  METHOD 

mlF°r  6aCh  Channel  to  be  filtered,  the  subroutine  erases 
locations  of  X,  and  multiplexes  the  filter  and  the 
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METHOD  (Contd.) 


data  channel  in  X,  starting  at  the  beginning'.  Note  that 
as  far  as  COOL  is  concerned,  X  is  a  complex  array  with 
data  in  the  real  part  and  filter  in  the  imaginary  part. 
COOL  is  called,  and  the  logic  of  COOLER  (q.v.)  is  used  to 
form  the  Fourier  transform  of  the  filtered  channel  in  X . 
COOL  is  called  again  to  get  back  to  the  time  domain,  and 
the  filtered  channel  is  written  on  the  output  tape. 

The  subset  label  is  copied  from  the  input  tape  to  the 
output  tape  at  the  beginning  of  the  subroutine. 


SEISMIC  DATA  LABORATORY 
ALEXANDRIA,  VIRGINIA 

DIGITAL  COMPUTING  SECTION 

A.  IDENTIFICATION 

Title :  Hilbert  transform  of  periodic  data 

COOP  Identification:  UES  G619  COOLHLBR 
Category:  G6  Time  Series  Analysis 

Programer :  E.  A.  Flinn  and  J.  F.  Claerbout 

Date :  23  September  1966 

B.  PURPOSE 

To  compute  the  Hilbert  transform  (quadrature  function) 
of  a  time  series.  Since  COOL  is  used,  the  time  series  is 
assumed  to  be  periodic  outside  the  range  of  definition. 

C  .  USAGE 

1.  Operational  Procedure:  This  is  a  FORTRAN  oj  subroutine, 
with  calling  sequence:  CALL  COOLHLBR(N,X) .  This  tub- 
routine  calls  COOL. 

2.  Parameters:  N  is  the  log  (base  2)  of  the  number  of 

data  points.  X  is  the  data,  dimensioned  at  least 
N 

2  in  the  calling  program,  and  typ^  complex  there. 

On  input,  the  real  da>ta  series  must  be 
stored  in  the  real  part  of  X,  and  the  imaginary  part 
must  be  zero . 

On  return,  the  real  data  series  is  stored 

N-l 

in  the  real  part  of  scaled  up  by  2  .  The  Hilbert 

transform  is  stored  in  the  imaginary  part  of  X,  also 

.  _  .  N-l 

scaled  up  by  2 
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3.  Space  Required;  Very  little  in  addition  to  the  array 

N+l 

for  data  whLch  requires  2  locations  m  the  call  mg 
program . 

4.  Temporary  Storage  Required:  None 


5. 

Printout:  None 

6. 

Error 

Printouts:  None 

7  . 

Error 

Stops:  None 

8. 

Input 

and  Output  Tape  Mountings:  Not  Applicable 

9. 

Input 

and  Output  Formats. 

Not  Applicable 

10. 

Selective  Jumps  and  Stops: 

None 

11. 

Timing 

:  Dominated  by  two 

calls  to  COOL 

12.  Accuracy:  The  data  is  returned  correct  to  ten  decimal 
places . 

13.  Cautions  to  User:  The  data  must  be  arranged  as  under 
(2)  above. 

Notice  that  as  far  as  this  subroutine 
is  concerned,  the  data  is  periodic  outside  the  range  of 
definition.  End  effects  may  cause  answers  which  the 
user  does  not  expect.  For  example,  if  the  input  is  a 
pure  sine  wave,  the  user  expects  the  quadrature  to  be 
a  pure  cosine.  Using  this  subroutine,  this  turns  out 
to  be  the  case  only  if  the  data  series  contains  an 
integral  number  of  cycles. 

14.  References  Writeup  of  UES  G612  COOL. 

D,  METHOD 

The  Hilbert  transform  of  a  function  has  a  Fourier 

L 

transform  which  is  (-1)  times  the  Fourier  transform  of 


R-2  0 


£>'  METHOD  jContd  ) 

th£  functior-  COtL  returns  the  real  and  imaginary  parts  of 
the  Founer  transform  of  a  function  calculated  from  zero 
to  2rr  so  that  the  real  part  is  symmetric  about  the  middle 
and  the  imaginary  part  is  antisymmetric. 

If  the  Fourier  transform  of  the  function  is  A+iB,  the 
Fourier  transform  of  the  Hilbert  transform  is  -B+iA..  All 
COOLHLBR  does  is  erase  the  second  half  of  the  Fourier 
transform  (the  part  from  rr  to  2rH  .  half-weight  the  end 

points,  and  call  COOL  again  to  transform  back  to  the  time 
domain . 

The  scale  factor  2*  1  comes  from  the  fact  that  COOL 


gives  the  unnormalized  transform. 
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SEISMIC  DATA  LABORATORY 
ALEXANDRIA ,  VIRGINIA 

DIGITAL  COMFUTLNG  SECTION 


A  IDENTIFICATION 

Title:  Fast  convolution  of  two  time  series  using  COOL 

COOP  Identification.  UE3  COOLVOLV 

Category  Time  series  analysis 

Programer  E.  A.  Flinn  and  D.  W.  McCowan 

Date:  23  September  1966 

B.  PURPOSE 

To  form  the  convolution  of  two  time  series,  not  by 
the  usual  polynomial  multiplication  algorithm,  but  by 
forming  the  two  Fourier  transforms  (using  COOL),  multi¬ 
plying  them  together,  and  transforming  back  to  the  time 
domain.  This  is  faster  than  the  usual  procedure  when 

LX'LF  »  4 ( 2N  +  1)  (LX  +  LF ) 

where  LX  is  the  data  series  length  LF  in  the  filter 
impulse  response  length  and  N  is  the  log  (base  2)  of 
LX  +  LF. 

C  .  USAGE 

1.  Operational  Procedure :  This  is  a  FORTRAN-63  sub¬ 
routine,  with  calling  sequence: 

CALL  COOLVOLV (LX  X , LF , F ) 

2 .  Parameters : 

X  is  the  data  series  to  be  convolved,  dimensioned  at 
least  2J+1  in  the  calling  program,  where  2  is  the 
smallest  power  of  two  larger  than  LX  +  LF  . 

LX  is  the  length  of  the  data  series  to  be  convolved. 
F  is  the  filter  to  be  convolved  with  X. 

LF  is  the  length  of  the  filter. 
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i  •  Space  Required:  300.^  Plus  arrays. 

4.  Temf-wf^ry  Locations  Required;  None  beyond  filling  out 
X  to  the  first  power  of  two  greater  than  LX  +  LF  . 

5  Alarms  or  Special  Printout ,  None 

6*  Error. Returns:  If  LX  +  LF  >  213,  LF  is  replaced  by  -LF 
and  control  is  returned  to  the  calling  program. 

7.  Error  Stops:  None 

8  .  Tape  Mountings:  None 

9.  Formats:  None 

10.  Jump  and  Stop  Settings:  None 

11.  Timing:  Dominated  by  two  calls  to  COOL  for  LX  +  LF 
points  each  time. 

Accuracy:  Gives  the  same  results  as  polynomial  multi¬ 
plication  to  ten  decimal  places. 

13.  Cautions:  None 

14 .  Configuration:  Standard  COOP 

15 •  References:  Writeups  of  COOL,  COOLCON.  and  COOLER 
D .  METHOD 

The  same  method  is  used  as  used  in  COOLCON. 
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APPENDIX  C  -  PROCEDURES 


FINITE  FOURIER  TRANSFORM  THEORY  AND  ITS  APPLICATION  TO  THE 
COMPUTATION  OF  CONVOLUTIONS,  CORRELATIONS,  AND  SPECTRA 


PROCEDURE  for  CALCULATING  A  CRo.q S 


SPECTRUM  AND  A  CROSS-CORRET.ATT on 


Dimension  X(2*LX+2),  CX(LX+1),  Y(2*LX+2),  CY(LX+1) 
Equivalence  (X,CX),  (Y,CY) 

Type  Complex  CX,CY 
LX  =  2 **N 

1)  Erase  2*LX+2  points  in  both  X  and  Y. 

2)  Read  channel  1  into  X  and  channel  2  into  Y. 

3)  Call  COOLER (N+ 1 ,  X ) 

Call  COOLER (N+l ,  Y) 

4)  Go  through  the  LX+1  complex  points  and  overlay  CX  (or  CY) 
vith: 

CX(I)  =  [CONJG (CX  ;)* CY(I)]/LX 
that  is, 

RerCX(I)  J  =  (RerCX(l)  l*Re(CY(l)  ]+imrcx(I)  ■’*im[CY(l)  ])/LX 

ImfCX ( I )  ]  =  (ReTcXd)  ]*Im[CY(l)  l-ImTCXd)  l*Re[CY(l)  ])/LX 

The  cross-spectrum  between  channel  1  and  channel  2  (which 
is  the  complex  conjugate  of  the  cross-spectrum  between 
channel  2  and  channel  1)  is  now  in  CX,  LX+1  points  in 
length.  The  co-spectrum  is  in  the  real  part  of  CX  and  the 
quad-spectrum  is  in  the  imaginary  part  of  CX . 

5)  To  get  the  cross-correlation,  fill  in  the  other  LX-1  points 
in  CX  and  call  COOL: 

DO  1  I  =  1 , LX-1 

1  CX ( LX+ 1+ 1 )  =  CONJG (CX(LX-I+1 ) ) 

CALL  COOL (N+l, CX, -1.0) 

The  cross-correlation  is  in  the  real  part  of  CX,  purely 
real  and  2*LX  points  in  length. 

NOTE:  CX  must  be  dimensioned  2*LX  if  the  cross-correlation 
is  to  be  calculated. 
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PROCEDURE  FOR  CALCULATING  AN  AUTO-SPECTRUM  AND  AN  AUTO- 
CORRELAT  ION - - - _ - - - 

Dimension  X(2*LX+2),  CX(LX+1) 

Equivalence  (X,CX) 

Type  Complex  CX,  CONTG 
LX  =  2**N 

1)  Erase  2*LX+2  points  in  X;  the  extra  complex  point  is  needed 
by  COOLER  to  eturn  the  point  at  the  folding  frequency. 

2)  Read  the  data  channel  into  X ( 1 .)  through  X(LX). 

3)  Call  COOLER ( N+ 1 ,CX)  .  The  Fourier  transform  of  X  and  the 
necessary  zeros  on  the  end  of  the  data  is  now  stored  in 
CX,  LX+1  complex  points  long,  representing  frequencies 
between  DC  and  the  folding  frequency. 

4)  Go  through  the  LX+1  complex  points  in  CX ,  and: 

CX ( I )  =  fCONJG (CX ( I ) ) *CX ( I ) ]/LX 
that  is, 

2  2 

Re[CX(I)1  =  (Re[CX(I)1  +  ImTcxtl)]  )/LX 
ImTCXd)  ]  =  0.0 

The  auto-spectrum  is  the  real  part  of  CX ,  purely  real 
and  LX+1  points  in  length. 

3)  To  get  the  auto-correlation,  fill  in  the  other  LX-1 
complex  points  in  CX  as  required  by  COOL  for  inverse 
transforms,  and  call  COOL: 

DO  1  I  =  1 , LX-1 
1  CX (X+ 1+ I )  =  CX (LX-I+ 1 ) 

CALL  COOL (N+l ,CX , -1 .0) 

The  auto-correlation  is  in  the  real  part  of  CX,  purely  real 
and  2* LX  points  in  length. 

NOTE:  CX  must  be  dimensioned  2*LX  if  the  auto-correlation 

is  to  be  computed. 
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PROCEDURE  FOR  CALCULATING 


THE  CONVOLUTION  OF  TW n  SERIES 

Dimension  X(L+2),  CX(Wl),  F(L+2),  CF^L+l) 
Equivalence  (X,CX),  (F,CF) 

Type  Complex  CX,CF,CONJO 
L  =  2**N 


L  here  is  the  next  power  of  2  larger  than  LX+LF, 
length  of  the  data  and  the  filter. 


the  combined 


1)  Erase  L+2  points  in  X  and  F. 

2)  Read  the  date  into  XU)  throu,h  XttX)  end  the  filter  impulse 
response  into  F ( 1 )  through  F{LF). 

3)  Call  COOLER (N.CX) 

Call  COOLER (N,CF) 


4)  Go  through  the  *SL+1  complex  points  in  CX,  and: 

CX(I)  =  rcX(I)*CF{l)  l/LX 
that  is, 

Re[CX(l) 1  -  (HercX(l)]*RerCF(I)l-i,flfcX(I)]*imrCF(I)])/LX 
I-nrcx(i)]  =  (RerCx(I)]*imrcF(I)]+RerCF(l)l*imrCX(l)-j)/LX 

The  Fourier  transformof  X  convolved  with  F  is  now  in  CX. 


5) 


Fill  in  the  rest  of  the  points  in  CX  as 
and  transform  back.  Note  again  that  if 
lution  is  desired  instead  of  che  Fourier 
must  be  dimensioned  L. 


needed  by  COOL, 
the  actual  oonvo- 
transform,  CX 


DO  1  I  =  1,  'jjL-l 
1  CX(>iL+I+l) 


CALL  COOL (N,CX, -1,0) 

The  convolution  of  X  with  F  is  now  in  the  reel  part  of  CX 
purely  real,  and  LX+LF-1  points  in  length. 
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